Computing intersections with polygons

Clipping and intersection functions can be used to extract trajectory segments that are located within an area of interest polygon.

import pandas as pd
import geopandas as gpd
from geopandas import GeoDataFrame, read_file
from shapely.geometry import Point, LineString, Polygon
from datetime import datetime, timedelta
import movingpandas as mpd

import warnings
warnings.filterwarnings('ignore')

print(f'MovingPandas version {mpd.__version__}')

MovingPandas version 0.8.rc1

gdf = read_file('../data/geolife_small.gpkg')
traj_collection = mpd.TrajectoryCollection(gdf, 'trajectory_id', t='t')

Clipping a Trajectory

help(mpd.Trajectory.clip)

Help on function clip in module movingpandas.trajectory:

clip(self, polygon, point_based=False)
    Return trajectory segments clipped by the given polygon.
    
    By default, the trajectory's line representation is clipped by the polygon.
    If pointbased=True, the trajectory's point representation is used instead, leading to shorter
    segments.
    
    Parameters
    ----------
    polygon : shapely Polygon
        Polygon to clip with
    point_based : bool
        Clipping method
    
    Returns
    -------
    TrajectoryCollection
        Clipped trajectory segments

xmin, xmax, ymin, ymax = 116.365035,116.3702945,39.904675,39.907728
polygon = Polygon([(xmin,ymin), (xmin,ymax), (xmax,ymax), (xmax,ymin), (xmin,ymin)])
polygon_gdf = GeoDataFrame(pd.DataFrame([{'geometry':polygon, 'id':1}]), crs=31256)

my_traj = traj_collection.trajectories[2]
intersections = my_traj.clip(polygon)
    
print("Found {} intersections".format(len(intersections)))

Found 1 intersections

ax = my_traj.plot()
polygon_gdf.plot(ax=ax, color='lightgray')
intersections.plot(ax=ax, color='red', linewidth=5)

<AxesSubplot:>

Clipping a TrajectoryCollection

Alternatively, using TrajectoryCollection:

clipped = traj_collection.clip(polygon)
clipped

TrajectoryCollection with 2 trajectories

clipped.plot()

<AxesSubplot:>

Computing intersections

help(mpd.Trajectory.intersection)

Help on function intersection in module movingpandas.trajectory:

intersection(self, feature, point_based=False)
    Return the trajectory segments that intersects the given feature.
    
    Feature attributes are appended to the trajectory's dataframe.
    
    By default, the trajectory's line representation is clipped by the polygon.
    If pointbased=True, the trajectory's point representation is used instead, leading to shorter
    segments.
    
    Parameters
    ----------
    feature : shapely Feature
        Feature to intersect with
    point_based : bool
        Clipping method
    
    Returns
    -------
    TrajectoryCollection
        Ssegments intersecting with the feature

polygon_feature = {
    "geometry": polygon,
    "properties": {'field1': 'abc'}
}

my_traj = traj_collection.trajectories[2]
intersections = my_traj.intersection(polygon_feature)
intersections

TrajectoryCollection with 1 trajectories

intersections.to_point_gdf()

	id	sequence	trajectory_id	tracker	geometry	intersecting_field1
t
2009-02-04 10:43:04.222205	3225	773	3	2	POINT (116.36799 39.90467)	abc
2009-02-04 10:43:05.000000	3226	774	3	2	POINT (116.36798 39.90471)	abc
2009-02-04 10:43:07.000000	3227	775	3	2	POINT (116.36795 39.90480)	abc
2009-02-04 10:43:09.000000	3228	776	3	2	POINT (116.36793 39.90487)	abc
2009-02-04 10:43:11.000000	3229	777	3	2	POINT (116.36794 39.90495)	abc
...	...	...	...	...	...	...
2009-02-04 10:48:14.000000	3390	938	3	2	POINT (116.36535 39.90589)	abc
2009-02-04 10:48:15.000000	3391	939	3	2	POINT (116.36525 39.90589)	abc
2009-02-04 10:48:16.000000	3392	940	3	2	POINT (116.36516 39.90589)	abc
2009-02-04 10:48:17.000000	3393	941	3	2	POINT (116.36507 39.90589)	abc
2009-02-04 10:48:17.399995	3393	941	3	2	POINT (116.36504 39.90589)	abc

170 rows × 6 columns